/*********************************/
* Table 3 results
* number of procedures and per-capita utilization rates (per 100,000 per year)
* for income quintile
* inputs:
*	indata = the cohort tables saved in the output folder
*	population_nsw = the population table (abs)
*	population_on = ontario population table (for standardization);
/**********************************/

/* convert populations */
* nsw population by income is in codes.population_nsw; 
data codes.population_nsw; 
	set codes.population_nsw;
	format sa2_name $45.; 
run; 

* ontario values; 
proc import out = codes.population_on_income datafile = 'C:\Users\kcha0642\Documents\NSW-ON-NY\on-ny-nsw-kcha0642\on-ny-nsw\population_on_income' 
	dbms = xlsx replace; 
	getnames = yes; 
run; 
data codes.population_on_income; 
	set codes.population_on_income; 
	format agegp $7.; 
run; 
data codes.population_on_income; 
	set codes.population_on_income; 
	sex_c = put(sex,$1.); 
run; 
data codes.population_on_income; 
	set codes.population_on_income; 
	format sex_c $HIRDSEX_1044_33.; 
run; 

%macro table_three (indata, population_nsw, population_on, label);

**********************************;
/* sa2 level (2016) and income quintiles */ 
proc sql; 
	create table temp as
	select a.SA2_2011_CODE, age_group, sex, count(distinct a.recnum) as n_levels 
	from &indata as a
	group by SA2_2011_CODE, age_group, sex;

* add abs correspondance file ; 
proc sql; 
	create table temp_cohort as
	select a.n_levels, a.SA2_2011_CODE as sa2_2011, b.sa2_maincode_2016 as sa2_2016, age_group, sex,
			b.sa2_name_2016 as sa2_2016_name, b.ratio
	from temp as a
	left join codes.sa2_cor as b
	on a.SA2_2011_CODE = b.sa2_maincode_2011; 

data temp_cohort; 
	set temp_cohort; 
	converted = n_levels * ratio; 
run; 

proc sql; 
	create table temp as
	select a.sa2_2016, a.sa2_2016_name, age_group, sex, sum(converted) as n_procedures
	from temp_cohort as a
	group by sa2_2016,  a.sa2_2016_name, age_group, sex; 

* add income labels; 
proc sql; 
	create table temp2 as
	select a.*, b.quintile_personal
	from temp as a
	left join codes.income_sa2 as b
	on a.sa2_2016 = b.sa2; 

**********************************;
* number of procedures; 

title "Number of procedures per income quintile"; 
proc tabulate data = temp2;
	class quintile_personal;
	var n_procedures; 
	table quintile_personal ALL, n_procedures*(sum pctsum);
run;

**********************************;
* age and sex adjusted; 

proc sql; 
	create table temp3 as
	select age_group, sex, quintile_personal, sum(n_procedures) as nsw_procedures
	from temp2
	group by age_group, sex, quintile_personal;

* procedures per year ; 
*** 5 years 3 months = 63 months; 
data temp3 ;
	set temp3;
	rate_nsw = nsw_procedures / (63/12); 
run; 

* nsw-population specific rates; 
proc sql; 
	create table nsw_incomes as
	select age, sex, quintile_personal, sum(count) as nsw_pop
	from codes.population_nsw
	group by age, sex, quintile_personal; 

proc sql; 
	create table temp4 as
	select a.*, b.nsw_pop
	from temp3 as a
	left join nsw_incomes as b
	on a.age_group = b.age and a.sex = b.sex and a.quintile_personal = b.quintile_personal; 

data temp4; 
	set temp4; 
	age_sex_specific_rate = rate_nsw / nsw_pop; 
run; 

* count in standard (ontario) population; 
proc sql; 
	create table temp5 as
	select a.*, b.sum_pop as ontario_population
	from temp4 as a
	left join &population_on as b
	on a.sex = b.sex_c and a.age_group = b.agegp and a.quintile_personal = b.incquint; 

* age and sex population rates; 
data temp5; 
	set temp5; 
	standard_pop_rate = age_sex_specific_rate*ontario_population; 
run; 

* export these tables for country comparison; 
proc export 
  data = temp5 
  dbms=xlsx 
  outfile="Z:\LowValueCare\APDC\on-ny-nsw\&label. income adjusted counts" 
  replace;
run;

* total utilisation rates; 
proc tabulate data = temp5 out = total_counts; 
	class quintile_personal; 
	var standard_pop_rate ontario_population; 
	table quintile_personal ALL, standard_pop_rate ontario_population; 
run; 

data total_counts; 
	set total_counts; 
	asr = (standard_pop_rate_Sum / ontario_population_Sum)*100000;
run; 

title "Total age-sex adjusted rates by income quintile"; 
proc print data = total_counts(keep=quintile_personal asr); 
run; 

%mend table_three; 




